Collective beating of artificial microcilia 
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We combine technical, experimental and theoretical efforts to investigate the collective dynamics 
of artificial microcilia in a viscous fluid. We take advantage of soft-lithography and colloidal self- 
assembly to devise microcapets made of hundreds of slender magnetic rods. This novel experimental 
setup is used to investigate the dynamics of extended cilia arrays driven by a precessing magnetic 
field. Whereas the dynamics of an isolated cilium is a rigid body rotation, collective beating results 
in a symmetry breaking of the precession patterns. The trajectories of the cilia are anisotropic 
and experience a significant structural evolution as the actuation frequency increases. We present a 
minimal model to account for our experimental findings and demonstrate how the global geometry 
of the array imposes the shape of the trajectories via long range hydrodynamic interactions. 
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Beating cilia are frequently encountered in nature 
to achieve propulsion, and to pump fluid at micro- 
scopic scales. Prominent examples include microor- 
ganisms, such as paramecia [lj, algae, such as Vol vox 
colonies [2j [3], and ciliated epithelial tissues, which di- 
rect flow of mucosa and fluids over macroscopic scales, 
see e.g. [4]- Recently, these propulsion mechanisms have 
sparked much interest in two distinct scientific communi- 
ties. From a technological perspective, the past few years 
have witnessed the quick development of cilia-inspired 
microactuators, aimed at transporting and mixing fluids 
in microfluidic channels. These works focused mainly on 
the flow profile induced by the asymmetric actuation of 
the artificial cilia [5 , 6j. From a more fundamental per- 
spective, the synchronization between beating cilia has 
motivated a surge in theoretical works [7J, echoed by a 
few experimental studies on artificial |8J and biological 
setups |9J. So far, these experiments have been dedi- 
cated to two-body problems. One of the main motiva- 
tions for these works was to understand the role of the 
hydrodynamic coupling in the coordination of cilia beat- 
ing in a viscous fluid. Most of the experiments and mod- 
els have focused on the phase dynamics of these coupled 
active/actuated systems. Conversely, we report on the 
beating amplitude of actuated microcilia carpets. 

In this letter, we investigate the dynamics of magnetic 
microrods driven by an external precessing field in a vis- 
cous fluid, see Fig. [I] We show that the long-range hy- 
drodynamic interactions result in an unexpected symme- 
try breaking of the beating trajectories, Fig. [2]A. To ad- 
dress this many-body problem, we propose a novel design 
strategy for the fabrication of extended magnetic micro- 
carpets, which we briefly describe below. In addition, we 
introduce a minimal but quantitative model to rational- 
ize our experimental findings. It indicates that the shape 
and orientation of the trajectories reflect the large-scale 
geometry of the cilia array. 

Magnetic microcilia carpets and experimental set-up. 
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FIG. 1: A: Snapshot of a square array of self-assembled col- 
loidal filaments driven by a precessing magnetic field. The 
pitch of the array is p = 30 /mi. The filaments are organized 
into a rectangular carpet inside a microchannel. B, C, D and 
E: The filaments are formed by self-assembly of superpara- 
magnetic colloids. A fraction of them aligns with magnetic 
anchoring sites at the bottom of the chamber. The exceden- 
tary filaments are rinsed away, leaving an array of filaments 
organized in the geometry imposed by the magnetic template. 



Over the last four years, numerous fabrication strate- 
gies have been proposed to pattern microchannel surfaces 
with field-responsive magnetic microcilia arrays [5j|6j[T0]. 
However, the accurate control of extended carpet geome- 
tries and the actuation via spatially homogeneous fields 
have not been achieved simultaneously. We combined 
soft lithography techniques and colloidal self-assembly to 
overcome these technical obstacles. By doing so, we man- 
aged to make prototypal magnetic cilia arrays organized 
into tunable arrangements over large scales, Fig.[T]A. The 
fabrication of the microfilaments themselves, via self- 
assembly of superparamagnetic colloids, has been previ- 
ously described [llj. However, the organization of these 
structures into controlled geometries remained virtually 
impossible. To achieve spatial ordering, we guide the self- 
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assembly of colloidal chains, using a soft-lithographied 
magnetic template sketched in Fig. [TJ3. In brief, using 
conventional replica molding techniques, we pattern a 
PDMS elastomer sheet with square lattices of square 
holes (width: 6 /im, depth: 3 jam and pitch: 30 /im, or 
40 /im). In each hole, a single paramagnetic colloid (Dyn- 
abeads, diameter: 3.5 /im) is deposited via capillary as- 
sembly following |12| . The resulting template, Fig. 
is used to seal a PDMS microchannel (height: 150 /im, 
width: 400 /im, length 1 cm). The channel is filled with 
colloidal particles of radius a = 375 nm and magnetic sus- 
ceptibility x ~ 1 (Ademtech). The colloids are diluted 
to a volume fraction ~ 0, 125% in an aqueous solu- 
tion containing 0.1 wt% polyacrylic acid (Mw=250 000, 
Sigma) and 0.1 wt% nonyl phenol ethoxylate (surfactant 
NP10, 0.1 wt%, Sigma), Fig.[l]3. The flow is then stopped 
with pneumatic PDMS valves [13J, and a 24 mT vertical 
magnetic field is applied immediately after the filling of 
the channel. The magnetic dipolar interactions between 
the paramagnetic particles cause them to organize into 
single stranded chains, aligned with the direction of the 
field and spread all over the channel, Fig. [Tp. The poly- 
acrylic acid molecules adsorbed on their surface link the 
colloids irreversibly jTT] . The excendentary filaments are 
subsequently rinsed away with a 50 wt% water-glycerol 
mixture (viscosity: rj = 5mPa.s), Fig. [l^. The resulting 
microcilia carpet is actuated by a set of three perpendic- 
ular Helmoltz coils, generating a 3D magnetic field, ho- 
mogeneous over the whole sample [T4l H5j . The internal 
dipolar interactions along each filament tend to align its 
main direction with the magnetic field, see Fig. [If and 
supplementary movie. We filmed simultaneously ~ 50 
cilia at 30 frames per second on an inverted Nikon mi- 
croscope (TE2000) with a 20x objective, Fig [l]A. Inci- 
dentally,note that this prototyping technique allows in 
principle for a vast range of geometries, which can ex- 
tend up to several millimeters. 

Collective beating of magnetic cilia. We consider the 
simplest istotropic 3D actuation cycle: the magnetic field 
precesses around the vertical axis at a constant angular 
velocity cj, keeping a constant angle 0b = 15° with the 
z-axis, Fig. [Tp. We first recall that in this small incli- 
nation limit, an isolated cilia responds linearly to the 
driving field, and the position of its tip follows a circular 
trajectory at a constant angular speed uo [15J. This is 
not what is observed with the cilia carpets. All the cilia 
do move synchronously at a constant angular speed lj, 
but the tip trajectories break the rotational symmetry 
of the driving, as shown on the pictures in Fig. [2]A. The 
increase in the driving frequency results in the stretch- 
ing of the beating patterns. In addition, the angle, ip, 
between the main axis of these anisotropic trajectories 
and the x-axis decreases with oj (the x-axis is here de- 
fined by the channel orientation). Changing the sign 
of the actuation results in mirrored trajectories. In or- 
der to quantify these structural changes, we fitted the 
trajectories by ellipses, and plotted the variation of the 
anisotropy, and of the orientation with uj in Figs. [2)3 and 
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FIG. 2: A: Superimposed pictures revealing the trajectories 
followed by the tip of the filaments. Lattice pitch: 30 /im, 
field amplitude: B — 8mT. B, C, D and E: symbols: experi- 
mental data. Full line: theoretical prediction form the mean 
field model. Blue filled symbols (resp. red open symbols) cor- 
respond to lattices with a p = 40 /im pitch (resp. p = 30 /im). 
Field amplitudes: B — 4 mt and 8 mT. B: Anisotropy as a 
function of the rescaled angular velocity. Straight line: Best 
fit from the mean- field model (r x = 0.19r, r y = 0.43r). C: 
Orientation of the major axis as a function of the rescaled 
angular velocity. D: Variations of the major axis and minor 
axes with ojt, and comparison with the variation of the tra- 
jectory radius for an isolated cilium (black diamonds). Inset: 
Close up on the low speed regime. The major axis displays 
non-monotonic variations. E: Variation of the time-averaged 
phase lag between the tip and the field orientation and com- 
parison with the same phase lag for an isolated magnetic cil- 
ium (black diamonds). 



C. We defined the anisotropy as (7*1 — r*2)/ri, where r\ 
and T2 are respectively the major and the minor axis 
of the ellipses. The dispersion of the data mainly orig- 
inates from the polydispersity of the filament lengths: 
t = 135 jam ± 10 jam. The error bars correspond to the 
maximal deviation from the mean values. Before going 
further, let us introduce the orientational relaxation time, 
r, of an isolated rod. r is defined by the force balance be- 
tween the magnetic and viscous forces acting on the cilia 
upon a sudden change in the field orientation. These 
forces scale as: / mag ~ (axB) 2 / (/ao£) and f v ~ (±(^/ r ) 
respectively, where (± = Atttj / [\n(£ / a) + 1/2] is the rod 
friction coefficient in a slender body approximation[16j. 
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We measured r according to the procedure introduced 
in [15| , which yields r=(4.37 10~ 3 )/B 2 for our experi- 
ments, where the length of the filaments has been kept 
constant, r is typically of the order of a few seconds, ind 

First, we note that all the measurements obtained 
for different field amplitudes, B, collapse on the same 
master-curves when we rescale the angular speed uo by r, 
Figs. [2|3, C, D and E. This implies that cjt, the so-called 
Mason number, is the only dimensionless parameter in- 
volving the driving speed in this problem. At low speed, 
ujt <C 1, the anisotropy increases linearly with co>, and ijj 
slowly decays from a finite value close to 45°. At high 
speed, ujt ^> 1, the anisotropy plateaus to a value close 
to 0.5 and the trajectories align with the x— axis. In addi- 
tion, we compare the major and minor axes individually 
to the radius, ro, of the circular trajectory followed by 
an isolated cilium, see Fig. [2p. In fact, the variations of 
ro and V2 cannot be distinguished: they both start from 
£ sin Ob and decay to above a well-defined crossover at 
ujt = 1. Conversely, ri, remains close to the static value 
£ sin Ob up to ujt ~ 5, and then decays to with the 
same slope as r<i- Looking more closely at the low speed 
regime, we notice that r\ undergoes non-monotonic vari- 
ations and reaches its maximal amplitude at ujt ~ 1, 
see the inset in Fig. [2p. These observations reveal that 
the cilia have more than one relaxation time when beat- 
ing within a carpet. At this point, we can anticipate 
on our theoretical model and infer that those two relax- 
ation times arise from geometrical corrections to r, since 
all the data collapse with ujt. To close this overview, 
we point that the time-averaged phase lag, (</>), between 
the rod and the field orientation is significantly reduced 
compared to the case of a single cilium. Nonetheless, the 
shape of the curve {(j)){uor) is conserved, see Fig. [2^1. 

Theoretical model and physical interpretations. To ex- 
plain quantitatively the symmetry breaking of the trajec- 
tories, we introduce a minimal far-field model. The rect- 
angular cilia carpet is modeled by a 27V x 2M lattice of 
pointwise particles located at = p(i e x +j e y ) +r^- (t), 
at a distance £ from a solid wall, see Fig. [3] Each particle 
is driven by an external time-dependent force f(R^,t). 
Note that we disregard the effect of the magnetic cou- 
pling between the cilia. The velocity of the particle 
is related to the force acting on all the other particles via: 



dt^ij — G(R^- — R n>m ) • f(R n?m ,t), 



(i) 



where, for sake of simplicity, the hydrodynamic coupling 
between the particles is described by the Blake-Oseen 
tensor G. G is the Green function of the Stokes equation 
associated to a force monopole oriented parallel to a solid 
wall in a viscous fluid [17). By definition, G(0) = C _1 2, 
is the isotropic mobility tensor for an isolated particle. 
We now perform a mean field approximation and as- 
sume that all the cilia follow identical trajectories in a 
synchronous manner: r^-(t) = r(t). This approxima- 
tion is justified by our experiments: we did not observe 
any spatial heterogeneities in the phase of the cilia tips. 



Assuming that £ > p, at leading order in p/£, we have 
G(R) = ^ B BIt, for R ^ 0. After some algebra, 
Eq{T]then reduces to the equation of motion for a single 
anisotropic particle driven by a time dependent force: 



d t r = (aI + Pe x e x ) • f(r,t), 



(2) 



where the two effective mobility coefficients are given by: 
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To compute the trajectories, we now need to specify 
f(t). As we consider only small field inclinations, the 
magnetic actuation is well approximated by a rotating 
harmonic trap: f(t) = Vk [r — i*o(t)] 2 , with i*o(t) = 
£ sin Ob (cos ujte x -\- sin ojte y ) [15J. Note that relaxing 
this harmonic hypothesis for the force and the geometric 
condition i > p does not qualitatively change our predic- 
tions, as Eq. [2]retains the same form. In the present case 
Eq. [2] corresponds to the equation of motion for two un- 
coupled ID overdamped harmonic oscillators driven by 
sinusoidal forces. They are readily solved, and the beat- 
ing trajectories r(t) = [x(t),y(t)] are: 



x(t) 



?sin Ob 



y/1 + (UT x y 



: cos [uot + arctan^r^)] (5) 



, x IsinOB . r / m / _ \ 

y(t) = sin [out + arctan(o;Ty)] (6) 

y/l + (uJT y ) 2 

where r x = l/[k(a + 0)\ and r y = l/(ka) are the two 
relaxation times of the particles, along the x and y di- 
rections respectively. Defining the major and the mi- 
nor axes as the extrema of r(t), we can compute nu- 
merically ri(uj) for i = 1,2, i/j(uj) and A0(cj). To test 
our theoretical predictions, we first determined the two 
free parameters r x and r y by fitting the anisotropy data, 
Fig. [2^, and then used the same parameters to calcu- 
late the major/minor axes, the inclination and the phase 
lag. As shown in Fig. [2j our model yields excellent agree- 
ment with the experiments. This unambiguously proves 
that the anisotropic trajectories originate only from the 
hydrodynamic coupling between the cilia, and not from 
their magnetic interactions. 

Moreover, Eqs. [5] and [6] imply that the symme- 
try breaking of the trajectories reflects the large scale 
anisotropy of the cilia carpet. The trajectories are 
anisotropic if r x ^ r y , or equivalent ly if f3 ^ 0. This 
mobility coefficient vanishes for symmetric (square) car- 
pets only: indeed it arises from the interactions between 
particles located at a distance larger than the width of the 
cilia carpet, see Eq(4] This coupling enhances the par- 
ticle mobility along the ^-direction, thereby effectively 
stretching the beating trajectory. 
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FIG. 3: A: Sketch of geometry used in the theoretical model. 
B: Two synchronous rotators driven by a harmonic trap. C: 
Direction of the induced forces in the low-speed limit. The 
trajectories are effectively sheared by the hydrodynamic cou- 
pling. D: Same picture in the high-speed regime 



More quantitatively, the minor and the major axes are 
reached at £*, defined as d t r 2 (t)\ t * = 0. At low speed, 
expanding Eqs. [5] and [6] at first order in cjt, we obtain 
dtT 2 (t) = —2oj(r x /Ty)ujTCOs(2ujt) : which implies that the 
trajectories are oriented at ujt* = 45°. Interestingly, this 
orientation is independent of the two relaxation time val- 
ues, provided that r x ^ r y . We can provide more phys- 
ical insight into this purely geometric result by looking 
at the structure of Eq. [2] The equation has the same 
structure for all carpet sizes and aspect ratio: under- 
standing the two-body problem is sufficient to account for 
the structural change of the beating patterns. Eq. [2] can 
be thought of as the equation of motion for two hydro- 
dynamically coupled synchronous particles, see Fig. [3j3. 
When ujt <C 1, the particle closely follows the trap mo- 
tion, and the force acting on the cilia is tangent to the 
trajectory. Therefore, the force induced by the rotation 
of the second cilia on the first one, /3(f -e x )e x , is extremal 
for ujt = : |mod(7r), and vanishes for ujt = 0mod(7r). 
This results in a net shear of the initially circular tra- 
jectory, as depicted in Fig.[3p. Consequently, the trajec- 
tories are expected to be stretched at a 45° angle from 
the x-axis, which is precisely what we observe in our ex- 
periments, Fig. |2p. We emphasize that this last picture 
is generic, and does not depend on the specifics of the 
driving force. Since the amplitude of the induced force 



scales as rcj, the stretching of the trajectory is expected 
to grow linearly with the driving speed. Again this qual- 
itative prediction is in agreement with our measurements 
and our asymptotic analysis. Indeed, estimating x and 
y at t* = 7t/(4cj), we obtain r 2 = 1 — (—iy(r x /r y )(ujr), 
with i = 1,2. We thus correctly predict that the major 
axis increases with the driving frequency, while the minor 
axis decreases, Fig(2]D and inset. As a result the corre- 
sponding anisotropy {r\ —r2)lr\ = (r x /r y )(ujr) increases 
linearly with uj, in agreement with our experimental data 
shown in Fig. [2j3. 

In the opposite high speed limit, ujt ^> 1, the phase 
lag of the driven overdamped particles reaches tt/2: the 
positions at which the induced forces are extremal are ro- 
tated by a 7r/2 angle with respect to the low-speed value. 
Therefore, as sketched in Fig.[3p, the hydrodynamic cou- 
pling results in a net stretching of the trajectories along 
the x-axis. This simple picture is confirmed by Eqs (5] 
and [6] For ujt ^> 1, they reduce to the parametrization 
of an elliptic trajectory aligned with the x-axis. Both r\ 
and T2 decrease as 1/cj, since the particle is driven faster 
than it can respond to the field: r\ ~ t sin 6 b / \ujt x ) and 
V2 ~ £ sin 6 b I \ujt x ) . For ujt ^> 1, i/j decays to 0, and the 
anisotropy plateaus to a constant value given by the ra- 
tio between the two relaxation times, r x /r y . This again 
confirms our experimental findings, Fig (2^ and C. In ad- 
dition, ri increases for ujt <C 1, whereas it decreases for 
ujt ^> 1, reaching a maximal value close to ujt — 1 as 
shown in the inset in Figj5J 

In conclusion, we combined technical, experimental 
and theoretical efforts to investigate the collective dy- 
namics of actuated microcilia in a viscous fluid. By doing 
so, we uncovered unexpected morphologies in the beating 
trajectories. We showed that this counter- intuitive dy- 
namic response is a consequence of the long range hydro- 
dynamic interactions. Importantly, we found the shape 
of the precession patterns to be chiefly selected by the 
large scale geometry of the carpets. 

We thank Sandrine Ngo for fruitful interactions, and 
Eric Lauga for a critical reading of the manuscript. 
We acknowledge support by C'Nano IdF, Sesame lie de 
France and Paris emergence. 



[1] M. Sleigh, The Biology of Cilia and Flagella (Pergamon 
Press, 1962). 

[2] L. Kirk, Molecular- genetic Origins of Multicellularity and 
Cellular Differentiation (Cambridge University Press, 
1998). 

[3] K. Drescher, R. E. Goldstein, and I. Tuval, Proc. Nat. 

Acad. Sci. USA 107, 11171 (2010). 
[4] K. Sawamoto, H. Wichterle, O. Gonzalez- Perez, 

J. Cholfin, M. Yamada, N. Spassky, N. Murcia, J. Garcia- 

Verdugo, O. Marin, J. Rubenstein, et al., Science 311, 

629 (2006). 

[5] M. Vilfan, A. Potocnik, and B. Kavcic, Proc. Nat. Acad. 
Sci. USA 107, 1844 (2010). 



[6] A. R. Shields, B. L. Fiser, B. A. Evans, M. R. Falvo, 

S. Washburn, and R. Superfine, Proc. Nat. Acad. Sci. 

USA 107, 15670 (2010). 
[7] R. Golestanian, J. Yeomans, and N. Uchida, Soft Matter, 

DOI: 10.1039/c0sm01121e (2011). 
[8] B. Qian, H. Jiang, D. A. Gagnon, K. S. Breuer, and T. R. 

Powers, Phys. Rev. E 80, 061919 (2009). 
[9] M. Polin, I. Tuval, K. Drescher, J. P. Gollub, and R. E. 

Goldstein, Science 325, 487 (2009). 
[10] J. den Toonder, F. Bos, D. Broer, L. Filippini, M. Gillies, 

J. de Goede, T. Mol, M. Reijme, W. Talen, H. Wilder- 

beek, et al., Lab on a Chip 8, 533 (2008). 
[11] C. Goubault, P. Jop, M. Fermigier, J. Baudry, 



5 



E. Bertrand, and J. Bibette, Physical review letters 91, 
260802 (2003). 

[12] L. Malaquin, T. Kraus, H. Schmid, E. Delamarche, and 

H. Wolf, Langmuir 23, 11513 (2007). 
[13] J.-C. Galas, D. Bartolo, and V. Studer, New Journal of 

Physics 11, 075027 (2009). 
[14] A. Babataheri, M. Roper, M. Fermigier, and O. du Roure, 

J. Fluid. Mech., in press (2011). 



[15] N. Coq, S. Ngo, O. du Roure, M. Fermigier, and D. Bar- 
tolo, Phys. Rev. E 82, 041503 (2010). 

[16] C. Brennen and H. winet, Annual Review of Fluid Me- 
chanics 9, 339 (1977). 

[17] J. Blake, Mathematical Proceedings of the Cambridge 
Philosophical Society 70, 303 (1971). 



